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Abstract 

Inspired by the quantum computing algorithms for Linear Algebra problems [HHL09, TS13] 
we study how the simulation on a classical computer of fhis type of "Phase Esfimation algo¬ 
rithms" performs when we apply if fo solve fhe Eigen-Problem of Hermitian mafrices. The 
resulf is a complefely new, efficienf and sfable, parallel algorifhm fo compufe an approxi- 
mafe specfral decomposition of any Hermitian mafrix. The algorifhm can be implemenfed 
by Boolean circuifs in O(log^n) parallel time wifh a fofal cosf of 0(n“+^) Boolean operations. 
This Boolean complexify mafches fhe besf known rigorous O(log^n) parallel time algorifhms, 
buf unlike fhose algorifhms our algorifhm is (logarifhmically) sfable, so furfher improvemenfs 
may lead fo practical implemenfafions. 

All previous efficient and rigorous approaches to solve the Eigen-Problem use randomiza¬ 
tion to avoid bad condition as we do too. Our algorithm makes further use of randomizafion in 
a complefely new way, taking random powers of a unifary mafrix fo randomize fhe phases of 
ifs eigenvalues. Proving fhaf a tiny Gaussian perturbation and a random pol 5 momial power are 
sufficient to ensure almost pairwise independence of the phases (mod 27r) is the main technical 
contribution of this work. This randomization enables us, given a Hermitian matrix with well 
separated eigenvalues, to sample a random eigenvalue and produce an approximate eigen¬ 
vector in O(log^n) parallel time and 0(n‘") Boolean complexity. We conjecture that further 
improvements of our method can provide a stable solution to the full approximate spectral 
decomposition problem with complexity similar to the complexity (up to a logarithmic factor) 
of sampling a single eigenvector. 


1 Introduction 

1.1 General 

The eigen-problem of Hermitian matrices is the problem of computing the eigenvalues and 
eigenvectors of a Hermitian matrix. This problem is ubiquitous in computer science and engi¬ 
neering, and because of its relatively high computational complexity imposes a high computa¬ 
tional load on most modern information processing systems. 

The Abel-Ruffini theorem implies there is no deterministic expression for the eigenvalues 
of a general matrix A, and for this reason eigen-problem algorithms must be iterative. This 
gives rise to a host of problems: these algorithms are hard to analyze rigorously, and often turn 
out to be unstable, and hard to parallelize. 

’School of Computer Science and Engineering, The Hebrew University, Jerusalem, Israel 
^Center for Theoretical Physics, MIT, Email: Ieldar@mit.edu 


1 



As with many other problems in computer science, one typically considers an approximate 
spectral decomposition of a matrix. Thus, given a matrix A, we are usually interested not in 
its exact eigenvalues / eigenvectors, which may be very hard to compute, (and possibly very 
long to describe once computed), but rather in an approximate decomposition: 

Definition 1. Approximate Spectral Decomposition - ASD(A, 5) 

Let A be some n x n Hermitian matrix. An approximate spectral decomposition of A, with accuracy 
parameter S = l/poly(n) is a set of unit vectors ||wi|| = 1 such that each Vi has \\Avi — 

\vi \\2 < 6, for some number \i e M, and 


y] Kvivf - A 




where |j Aj] is the operator norm of X. 

By standard arguments the above can be generalized to an arbitrary n x n matrix A, by 
considering the Hermitian matrix A^ A, in which case ASD(A^ A, 6) is an approximation of fhe 
singular vectors (and singular values) of A. We nofe fhaf fhe definifion of ASD fhen corresponds 
fo a "smoofh analysis" of mafrices: namely given inpuf A, we do nof find a specfral decom- 
posifion of A, buf rafher fhe decomposifion of a mafrix A', such fhat ||A — A'|| < S. We also 
poinf ouf, fhaf fhe definifion of ASD holds jusf as well in the case of nearly degenerafe matri¬ 
ces: we do not require a one-to-one correspondence with the eigenvectors of A, which can be 
exfremely hard fo achieve, buf rafher fo find some sef of approximafe eigenvectors, such thaf 
fhe corresponding weighfed sum of rank-1 projecfions form an approximafion of A. 

When one considers an algorifhm A for fhe ASD problem, one can examine ifs arithmetic 
complexify or boolean complexify. The arifhmefic complexify is fhe minimal size arifhmefic cir- 
cuif C (namely each node computes addifion/mulfiplicafion/division fo unbounded accuracy) 
fhaf implemenfs A, whereas fhe boolean complexify counfs fhe number of boolean AND / OR 
gafes of fan-in 2 required fo implement A. 

Given the definition above, and following Demmel ef al. [DDH07] we consider an algo¬ 
rithm A to be log-stable (or stable for shorf), if fhere exisfs a circuif C fhaf implemenfs A on 
n X n matrices, and a number t = 0(log(n)), such that each arithmetic computation in C uses 
at most t bits of precision, and outpuf of fhe circuif deviates from fhe oufpuf of fhe arifhmefic 
circuif by af mosf 1/poly(n). We nofe fhaf when an algorifhm is stable fhen ifs boolean com¬ 
plexify is equal fo ifs arifhmefic complexify up fo a facfor 0(log(n)). If, however, an algorifhm 
is unstable fhen ifs boolean complexify could be larger by a facfor of up fo n. In fhe sfudy of 
practical numerical linear algebra algorifhms, one usually identifies algorifhms fhaf are sfable 
wifh "pracfical", and algorifhms fhaf are nof sfable fo be impracfical. This usually, because 
fhe computing machines are resfricfed fo representing numbers wifh a number of bifs fhaf is a 
small fraction of fhe size of fhe inpuf. 

In terms of parallelism, we will refer fo fhe complexify class (see Definition 10) which 

is the set of all compufafional problems fhat can be solved by uniform Boolean circuifs of size 
poly(n) in time 0(log^(n)). Often, we will refer fo fhe class RNC*^*'^, in which fhe parallel 
circuif is also allowed fo accepf uniform random bifs. One would like an ASD algorifhm fo 
have minimal arifhmefic / boolean complexify, and minimal parallel time. Ideally, one would 
also like fhis algorifhm fo be sfable. 

1.2 Main Contribution 

Inspired by recent quantum computing algorithms [HHL09, TS13], we introduce a new per¬ 
spective on the problem of compufing fhe ASD fhaf is based on low-discrepancy sequences. 
Roughly speaking, low-discrepancy sequences are deferminisfic sequences which appear fo be 
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random, because they "visit" each small sub-cube the same number of times that a completely 
random sequence would up to a small additive error. 


Definition 2. Multi-dimensional Discrepancy 

For integer s, put P = [0,1)®. Given a sequence x = {xn)n=i, with Xn € P the discrepancy Dm{.x) 
is defined as: 


Dm{x) = sup 
BeB 


1 

-V0l{B) 


where XB{xn) is an indicator function which is 1 if Xn € B is the set of all s-products of intervals 
Y\i=i[ui,Vi],with [ui,Vi]{modl) C [0,1). 


We recast the ASD problem as a question about the discrepancy of a cerfain sequence relafed 
fo fhe inpuf mafrix. Specifically, given a Hermifian mafrix A wifh, say n unique eigenvalues 
{Ai}ig[n] the central object of inferesf is fhe sequence comprised of n-dimensional vecfors of 
eigenvalue residuals: 


S{A) = ({Ai • 1},..., {A„ • 1}), ({Ai • 2},..., {A„ • 2}),... ({Ai • M},..., {A„ • M}), 


where {a;} is fhe fracfional parf of a; € K, and M = poly(n) is some large infeger. S'(A) is 
hence a sequence of lengfh M in [0,1)". We would like S{A) fo have as small discrepancy as 
possible. Hence, in sharp confrasf fo previous algorifhms, insfead of fhe compufafional effort 
being concentrated on revealing "structure" in the matrix, our algorithm is actually focused on 
producing random-behaving d 5 mamics. 

The main applicafion of our approach presenfed in this paper is a new stable and parallel 
and stable algorithm for computing fhe ASD of any Hermifian mafrix. 


Theorem 3. For any Hermitian matrix 0 ^ A ^ 0.9/, and S < n we have ASD(A, (5) G 
with total boolean complexity 0{rP^fi. The algorithm is log-stable. 


The boolean complexify of our algorifhm is 0(n‘^+^). If however, one is inferesfed in sampling 
a uniformly random eigenvector, it can be achieved in complexity 0{rP). ^ 

1.3 Overview of the Algorithm 

To compute the ASD of a given mafrix A, we firsf consider a similar problem of sampling 
uniformly an approximate eigenvecfor of A, where fhe eigenvalues of A are assumed fo be 
well-separated. Clearly, if one can sample from fhis disfribufion in RNC^, fhen by fhe coupon 
collector's bound concafenafing 0(nlog(n)) many parallel copies of fhis routine, one can sam¬ 
ple all eigenvectors quickly wifh high probabilify. To do fhis, we require a definifion of a 
Hermifian mafrix fhaf is (5-separafed: 

Definition 4. S-separated 

Let Abe an n x n PSD matrix. A has a complete set of real eigenvalues Ai > A 2 > ... > A„ > 0. We 
say that A is S-separated if Xj — Aj+i > 5 for all j < n, and Ai < 1 — 5. 

Next, we introduce the notion of a separafing infeger w.r.f. a sequence of real numbers: 

Definition 5. Separating Integer 

Let A = (Ai,..., A„) G [0,1)". For a > 4 define 

Bout = [1-1/{In),-1 + 1/{In)] and Bin{a) = [-l/{an),l/{an)], 

A positive integer m is said to separate the k-th element ofXw.r.t. Bin, Bout if it satisfies: 
signifies the infimum over all constants c such that one can multiply two matrices in time at most n'^. 


3 






• {mXk} G Bin{a) 

• yj ^ k {mXj} ^ Bout 

and finally define fhe nofion of a separafing infeger w.r.f. a (5-separafed mafrix. 

Definition 6. A positive integer m is said to separate k in a S-separated matrix A iv.r.t. Bin, Bout, if 
m separates the k-th element ofC{A) w.r.t. Bin, Bout- 

Following is a skefch of fhe main sampling roufine. For complefe defails see Section 5. 
The roufine accepfs a separafing infeger m of fhe f-fh eigenvalue of a (5-separafed mafrix A, a 
precision paramefer 8 and refums a 8 approximation of fhe f-fh eigenvecfor of A. 

Algorithm 7. Filter(A, to, i5) 

1. Compute parameters: 

p = 24n^|’ln(l/(5)],(^ = 8“^/{2pm). 

2. Sample random unit vector: 

Sample a standard complex Gaussian vector v, set wq = z;/||r'||. 

3. Approximate matrix exponent: 

Compute a C Taylor approximation of denoted by U. 

4. Raise to power: 

Compute t/"* by repeated squaring. 

5. Generate matrix polynomial: 

Compute B — ^ i+u'^ ^ by repeated squaring. 

6. Filter: 

Compute w = 

7. Decide: 

Set z = A ■ w, in = argmaxjg[„] IwiI and compute c = Zig/wi^. If 

IIA • w — c • w|| < 38^/n 

return w, and otherwise reject. 


In words - the algorithm samples a random vector and then multiplies it essentially by the 
matrix B = {{I + e®"^™)/2)P. After this "filtering" step, it evaluates whether or not the resulting 
vector is close to being an eigenvector of A, and keep if if if is. To undersfand fhe behavior of 
fhe algorifhm, if is insighfful fo consider fhe behavior in fhe eigenbasis of A. 

w = ^ aiWi, 


where is an orfhonormal basis for A corresponding fo eigenvalues {Ai}jg[„]. If {mXi}, 

i.e. - fhe fractional parf of mXi, is very close fo 0, and {toA^} is ^ 2 In n/p far from 0 for all j f i, 
fhen affer multiplication by B and normalization, all eigenvecfors wj for j f i are affenuafed 
by factor 1/rf relafive fo Wi, and hence fhe resulfing vecfor is 1/n close fo an eigenvector of A^. 

Hence, a sufficienf condition on fhe number to fhaf would imply fhaf w = FilterjA, to, 5) 
is an approximation of fhe f-fh eigenvecfor is fhe following properly: {TOAi} is very close fo 0, 
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and for all j ^ i {mXj} is bounded away from 0. This corresponds fo fhe facf fhaf m separates i 
in A, as assumed. 

So fo sample uniformly an approximafe eigenvecfor, we would like fo call Filter {A, m, S) 
for m ~ U[M] for M — poly(ri) such fhaf m separafes i where i ^ U[n]. The main observafion 
here, is fhaf fhis condifion is safisfied if fhe sequence of residuals of infeger mulfiples of fhe 
eigenvalues S'(A) defined above has fhe aforementioned low discrepancy properfy. 

Mosf of fhe work in fhis sfudy is devofed fo achieving fhis properfy. Compufafionally we 
achieve low-discrepancy of S{A) simply by additive Gaussian perfurbafion prior fo calling fhe 
sampling roufine. We show fhaf if we perfurb a mafrix using a Gaussian mafrix £ of variance 
l/poly(n), fhen S'(d. -|- f) has discrepancy which is l/poly(n). Showing fhis is non-frivial be¬ 
cause arbifrary vecfors of eigenvalues Ai,..., A„ do nof generafe low-discrepancy sequences 
in general ^ , and on fhe ofher hand we are also severely limifed in our abilify fo perfurb fhe 
eigenvalues wifhouf deviating foo much from fhe original mafrix. This is fhe subjecf of our 
main fechnical fheorem 32, which may be of independenf inferesf: 

Theorem, (sketch) Let Abe an n x n Hermitian matrix, and £ be a standard Gaussian matrix. For 
any a > 0,b > 0 there exists M = M {a, b) — poly(n) such that w.p. at least 1 — n~^ the sequence of 
residuals of eigenvalue multiples of A + • £ of length M has discrepancy at most n~^. 

Perfurbing fhe inpuf mafrix has fhe addifional benefif of making sure fhaf A has a exacfly 
n unique eigenvalues wifh high probabilify. This follows from a breakfhrough fheorem by 
Nguyen, Tao and Vu [NTV16] which has provided a resolufion of fhis long-sfanding open 
problem, which was considered unproven folklore until fhaf poinf. This fheorem allows us fo 
handle general Hermitian mafrices wifhouf exfra conditions on fhe condifioning number of A 
/ ifs eigenvalue spacing. 

1.4 Prior Art 

There are numerous algorifhms for computing fhe ASD of a mafrix, relying mosf prominenfly 
on fhe QR decomposifion [TB97]. For specific f 5 rpes of mafrices, like fridiagonal mafrices much 
fasfer algorifhms are known [ReiOS], buf here we consider fhe mosf general Hermitian case. 
We summarize fhe sfafe of fhe arf algorifhms for fhis problems in ferms of fheir complexify 
(boolean / arifhmefic, serial / parallel) and compare fhem fo our own: 



Arithmetic 

Gomplexity 

Boolean 

Gomplexity 

Parallel Time 

Log-Stable 

Gomments 

Csanky 

0(n‘^+i) 


log^(n) 

NO 


Demmel et 
al. [DDH07] 


0(n-)(*) 

N/A 

YES 

* Conjectured by us, by 
modifying the algorithm. 

Bini et al., 
Reif [BP92, 
ReiOS] 



0(log^(n)) 

NO 

Working with r2(n) bit In¬ 
tegers 

New 

0(n‘^+i) 

0(n‘^+i) 

log^(n) 

YES 



Gomparing our algorifhm fo fhe besf known algorifhms, if is more efficienf by a facfor 

of n compared wifh Gsanky's algorifhm [Koz92]. Nofably, our algorifhm is complefely disjoinf 

^Consider for example the sequence of values 1,1/2, • • ■, 1/n. Then multiplying these numbers individually by a 
random integer m and taking the residual would map to 0 all values 1/i for which i\m, and there is a logarithmic 
number of these on average. This sequence of eigenvalues is well-separated, and at least potentially could arise from a 
random matrix. We show, however, that this is not the typical case. 
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from Csanky's techniques - which rely on computing explicitly high powers of the input ma¬ 
trix, and computes the characteristic pol 5 momial of the matrix using the Newton identities on 
the traces of those powers. This is an inherently unstable algorithm as it finds the eigenval¬ 
ues by approximating the roots of the characteristic pol 5 momial and small perturbation to the 
coefficients of the pol 5 momial may lead to large deviations of the roots. 

The algorithms of Demmel et al., Bini et al. and Reif, rely on efficient implementation 
of variants of the QR algorithm. Our asymptotic bounds are worse then Demmel et al. in 
terms of total arithmetic/boolean complexity, though we conjecture that this is an artifact of 
our proof strategy, and not an inherent problem (see the section on open problems), and in 
fact, a variant of the algorithm could probably achieve a boolean complexity of 0(n“). We 
note that the QR algorithm is not known to be parallelizable in a stable way, and hence the fast 
parallel algorithms of Bini et al. and Reif are not stable and probably impractical. In fact the QR 
decomposition has been shown, for standard implementations like the Given's or Householder 
method, to be P-complete [LMM99] assuming the real-RAM model. Thus, it is unlikely to be 
stably-parallelizable unless P = NC. ^ 

Thus, to the best our knowledge, our algorithm is the first parallel algorithm for the ASD of 
general Hermitian matrices that is both parallel and stable. In particular it achieves the smallest 
bit-complexity of any algorithm to date. We conjecture that our approach may present 

a practical and parallel alternative to computing the ASD. We dwell on this point a bit more in 
Section 8. 

1.4.1 Comparison to the power method / QR algorithm 

An arguably natural benchmark by which to test the novelty of the proposed algorithm is the 
iterative power-method for computing the eigenvalues of a Hermitian matrix. In this method, 
one starts from some random vector bo, and at each iteration k sets: 

. _ Abk 

~ WAbkW 

This method can achieve polynomially good approximation of the top eigenvalue in time 
which is logarithmic in n, for A with a constant spectral-gap. 

Both the power method and our proposed scheme are similar in the sense that they attempt 
to extract the eigenvectors of the input matrix directly. Also, if two eigenvalues are e-close in 
magnitude, for some e > 0, then they require essentially the same exponent of A in the power 
method, and of in our scheme to distinguish between them. However, the similarity 
stops here. We maintain, that the power method is both conceptually different, and for gen¬ 
eral Hermitian matrices performs much worse, in terms of running time, compared with our 
proposed algorithm. 

Conceptually, in the power method, we seek to leverage the difference in magnitude between 
adjacent eigenvalues in order to extract the eigenvectors. Qn the other hand, in our proposed 
scheme we recast the problem on the unit sphere SAt, where we are interested in the spacing 
of the residuals of integer multiples of the eigenvalues. Worded differently, our setting exploits 
the additive group structure of the eigenvalues modulo 1, whereas the power method distin¬ 
guishes between them multiplicatively. In the additive group setting, the advantage is that 
we can consider the discrepancy of the sequence of residuals, and analyze how quickly these 
residuals mimic a completely independent random distribution. Furthermore, in the additive 
setting there is inherent symmetry between the eigenvalues, as no eigenvalue is more likely to 

^We point out that the algorithm of Reif [Rei05] achieves a QR factorization in parallel time 0(log^(n)) in the arith¬ 
metic model, thus showing that QR is indeed parallelizable, but it relies on computations modulo large integers and 
therefore not stable and not practical. 
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be sampled than another. This allows for a natural parallelization of fhe algorifhm fo exfracf 
simulfaneously approximation of all eigenvecfors. 

The well-known QR algorifhm for eigendecomposifion [GVL96] is fhe de-facfo sfandard for 
compufing the ASD, and is, in a sense, a parallel version of fhe power-mefhod. Thaf algorithm 
applies an iterated sequence of of QR decomposifions: Af each sfep k we compufe (where 
Ai = A - fhe inpuf mafrix) 

Ak — QkRkj 


and then set 


Ak+l — RkQk- 


The algorithm runs in time 0{rA), by applying several pre-processing steps [GVL96], and the 
fast variant of Demmel ef al. in time 0(n“). However, as stated above, the QR decomposition 
which is at the core of fhese mefhods is nof known fo be sfably parallel. 


1.5 Open Questions 

We oufline several open questions that may be interesting to research following fhis work: 

1. Is if possible fo affain a serial run-fime of 0{'n!^) for fhis algorifhm? We conjecfure fhaf 
fhis is possible based on numerical evidence for a varianf of fhis algorifhm, yef we do 
nof have a proof of fhis facf. While nof direcfly improving on fhe besf previously known 
serial run-fime ([DDH07]), if would reduce fhe overall work performed by our parallel 
RNC*-^^ algorifhm fo mafch fhe work done by sfafe-of-fhe-arf serial run-time algorithm. 

2. What other linear-algebra algorithms can be designed using our methods ? We would 
like these algorithms to improve on previous algorithms in either the stability, boolean 
complexity, parallel run-time, or all these parameters together. 

3. Gould one reduce the number of random bifs required by fhe algorifhm? Currenfly - we 
show fhaf using 0{n?) random bifs - i.e. applying additive Gaussian perfurbafion resulfs 
in a mafrix whose eigenvalues seed a low-discrepancy sequence. However, can one do 
away with only 0{n) random bits - by applying a tri-diagonal perturbation to the matrix? 

4. Is our algorithm practical ? What is the actual run-time of fhe algorifhm on mafrices 
of "reasonable" size, and does if compare wifh sfafe-of-fhe arf parallel algorifhms? Qur 
numerical evidence suggesfs fhaf in practice our algorifhm may run much fasfer fhan fhe 
analytical asymptotic bounds we provide here. 
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2 Preliminaries 


2.1 Notation 

A random variable x distributed according to distribution T) is denoted hy x ^ V. For a matrix 
X, ||X|| signifies the operator norm of X. For a sef S, [/[S'] is fhe uniform disfribufion on S. 
For infeger M > 0 fhe sef [M] is fhe sef of infegers {0,1,..., M — 1}. For real number x, {x} 
denofes fhe fracfional parf of x: {x} = x—\x\. For a Hermifian nxn mafrix A, wifh eigenvalues 
> A 2 > ... > A„ jC(A) = (Ai,..., A„) S K" denofes fhe vector of sorted eigenvalues 
of A. For a measurable subsef S C M" vol(S) denofes fhe volume of S. $ is fhe empfy sef. GUE 
is fhe global unifary ensemble of random mafrices. N, Z, C signify fhe nafural, infeger, and 
complex numbers, respecfively. For a mafrix A, A^ is fhe Flermifian conjugafe-franspose of A. 
For number n > 0 In n denofes fhe nafural logarifhm, and logn denofes fhe binary logarifhm. 

tr^) is fhe sfandard Gaussian measure wifh mean rj and variance U (n) is fhe sef of nxn 
unifary mafrices. 


2.2 Definitions 


2.2.1 Hermitian matrices 

We repeaf again fhe definifion of fhe ASD due fo ifs imporfance: 

Definition 8. Approximate Spectral Decomposition - ASD(A, (5) 

Let A be some nxn 6-separated Hermitian matrix. An approximate spectral decomposition of A, with 
accuracy parameter 6 = l/poly(n) is a set of unit vectors ||ni|| = 1 such that each Vi has 

\\Avi — \iVi \\2 < 5, for some number Ai G C, and 


Kv^vf - A 

iG[ra] 


<< 5 , 


where || A|| is the operator norm of X. 

By standard arguments the above can be generalized to an arbitrary nxn matrix A, by 
considering the Hermitian matrix A^ A, in which case ASD(A^A, S) is an approximation of 
fhe singular vectors (and singular values) of A. We nofe fhaf if {r’i}ig[„] is an ASD(A, <5) for some 
mafrix A, and B is a mafrix such fhaf || A — B\\ < 6, fhen by fhe friangle inequalify {vi]i^[n] is 
an ASD(A, ^2(5). 


2.2.2 Complexity 

Definition 9. Let uj denote the infimum over all t such that any two nxn matrices can be multiplied 
using a number of products at most n*. 

The current best upper-bound on w is 2.372 due to Williams [Will2]. 

Definition 10. Class NC 

The class is the set of problems computed by uniform boolean circuits, with a polynomial number 

of gates, and depth at most O(log^n). 

We will require the following known facf: 

Fact 11. [Akl90] There exists an algorithm for sorting n numbers in time log(n), using n processors. 

Definition 12. Class RNC 

The class is the set of problems that can be computed by uniform boolean circuits, with a poly¬ 

nomial number of gates, accepting a polynomial number of random bits, and depth at most O(log^n). 
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For simplicity, we shall assume in this work that RNC circuits are allowed to accept <-bit 
numbers, sampled from a truncated Gaussian distribution, and discretized to t-bits of preci¬ 
sion. 

2.2.3 Stable Computation 

Following Demmel et al. [DDH07] we define the notion of log-stability as one where truncating 
each binary arithmetic operation to 0(log(n)) bits of precision doesn't change the result by 
much: 

Definition 13. {t, 5)-stable randomized computation 

Let C denote a randomized arithmetic circuit, and V be its output distribution supported on M". Let D 
denote the discretization of C tot bits as follows: each infinite-precision arithmetic operation is followed 
by rounding to t bits. Let V denote the output distribution of D. C is said to be {t, 6)-stable if 

Mx 3y, 'D{x) = ’Dfy) and ||a; — y\\ < 5. 

Definition 14. Log-stable computation 

Let C be a randomized arithmetic circuit that accepts n input numbers. C is said to be log-stable if for 
any 6 = l/poly(n) it is {t, 5)-stable for some t = 0(log(l/(5)). 


3 Additive Perturbation 

Matrix perturbation is a well-developed theory [SS90, GVL96] examining the behavior of eigen¬ 
values and eigen-vectors under additive perturbation, usually much smaller compared to the 
norm of the original matrix. While general eigenvalue problems are usually unstable against 
perturbation, for Hermitian matrices the situation is much better: the Bauer-Fike theorem states 
that the perturbed eigenvalues can only deviate from the original eigenvalues by an amount 
corresponding to the relative strength of the perturbation. 

In particular, when the perturbed matrix A is ^-separated one can compute an explicit es¬ 
timate for the behavior of the perturbed eigenvalues. We use here a quantitative estimate by 
[SS90]: 

Fact 15. Rayleigh quotient for well-separated eigenvalues 

Let Abe a 6-separatednxn LLermitian matrix with eigenvalues \i > A 2 > ... > \n, and corresponding 
orthonormal basis het E be an additive perturbation of A satisfying \Eij\ < efor all i,j. Let 

\i denote the i-th eigenvalue ofA-\-£. There exists a constant c > 0 satisfying: 

Vz e [n] \i = X,-\- v^Ev^ + Q, ICil < ce'^/S- 

In fact, if the perturbation E is GUE a stronger characterization is readily available: 

Corollary 16. Let Abe a S-separated n x n Hermitian matrix with eigenvalues {Ai}ig[„], and corre¬ 
sponding orthonormal basis {wileH- ^ GUE. Then the eigenvalues {A'}ig[„] of the perturbed 
matrix A' = Ae ■ E are distributed as follows: 

Vz € [u] Ai = (1 — a) • p{\i,£^) a ■ V -\- C,i, 0<a< 2“", |^i| < 16cn • /5 

for some distribution V. 

Proof. By Fact 15 the eigenvalues A' behave as 

A- = Xi +v^Evi + C,i, ICil < cmax|£ij|^/(5. 
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The standard Gaussian satisfies: 


Pm (kl > 4v^) < 2-'”. 

Thus, by the union bound we have that \Sij\ < 4yTi for all i, j w.p. af leasf 1 — 2“". Hence, 
w.p. af leasf 1 — 2“" we have: 

Vi e [n] |Ci| < cmax|£i jp/J < 16cn ■ /S, 

Suppose fhaf fhis is fhe case, and lef £' denofe fhe GUE mafrix, conditioned on having bounded 
enfries. We can wrife: 

£' = {l-a)-£ + a-V, 

where V is some disfribufion onnxn mafrices and 0 < a < 2“". The disfribufion £ is invarianf 
under unifary conjugation, i.e.: 

VuSC/(n) u£u^ = £ 

fhen 

Mu G U{n) u£'u^ = {1 — a) ■ £ + a - V , 

where V = uVu^ is some disfribufion on n x n mafrices. Hence, up to statistical distance at 
most a we can assume w.l.o.g. that Vi = e^. This implies that 

= Ai + (1 — a) • £i^i + a • 

where |Ci| < 16cn • /5. Since £i^i = /x(0, e^) then 

A- = (1 - a) • p(Ai,e^) + a • P' + Cz, ICd < 16cn • /5. 


□ 

The stability of eigenvalues has also been generalized fo eigenvecfors and even general 
invarianf subspaces [SS90] in fhe following sense: if fhere is a clusfer of eigenvalues fhaf 
is "well-separafed" from all ofher eigenvalues, fhen fhe orfhogonal projection onto fhe sub¬ 
space spanned by fhe corresponding eigenvecfors is likewise sfable under additive perfurba- 
fion whose scale is negligible compared fo fhe separafion of fhe eigenvalue of fhis clusfer from 
fhe resf of fhe specfrum. In particular, if a mafrix is 5-separafed, fhen ifs eigenvecfors are indi¬ 
vidually sfable as follows: 

Fact 17. [SS90] Rayleigh quotient for well-separated invariant subspaces 

Let Abe a 6-separated n x n Hermitian matrix with eigenbasis {vi\i^[n], imd 

Ai=A + £, Ill’ll <e. 

There exists an orthonormal basis ^ constant c > 0 satisfying 

Mi G [n] \\vi — 'c'll < ce^/5. 

Our interest in additive perturbation, however, is not confined jusf fo "sfabilify" argumenfs. 
In facf, our main reason for using perfurbation is fo cause a scattering of fhe eigenvalues. The 
firsf step of our algorithm in fact applies additive perturbation to provide a minimal spacing 
between eigenvalues. Recently Nguyen et al. [NTV16] have provided the first proof fhaf ap¬ 
plying additive perfurbation fo any Hermifian mafrix using a so-called Wigner ensemble, an 
ensemble of random mafrices fhaf generalize GUE, in facf causes fhe eigenvalues of fhe per- 
furbed mafrix fo achieve a minimal inverse pol 5 momial separafion. We sfafe fheir resulf: 
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Lemma 18. [NTV16] Minimal eigenvalue spacing 

Let Mn = Fn + e ■ where is a real symmetric matrix, ||Fn|l 2 ^ 1/ e = rT^ for some constant 
7 > 0, and Xn is GUE. Let Ai > A 2 > ... > A„ denote the eigenvalues of Mn, and put ai = Xi — A^+i 
for all i < n. Then for any fixed ^ > 0 there exists B = _B(e) > 0, such that 

maxP (ai < n~^) = 0(n~^). 

iG[ra] 

In particular for any A > 0 there exists B > 0 such that 

P ( minai > | = 1 — 0(n~^). 

VeN J 

Using the lemma above we define: 

Definition 19. For any S = 1/poly (n), let B*(S) denote the largest number B > 0 such that for every 
Fn the matrix Mn = + SXn satisfies: 

P I minoi > n~^] > 0.99 

V*6H J 


4 Low-Discrepancy Sequences 


4.1 Basic Introduction 


Low discrepancy sequences (or "quasi-random" sequences) are a powerful fool in random 
sampling mefhods. Roughly speaking, fhese are deferminisfic sequences fhaf visif any "rea¬ 
sonable" subsef B a number of fimes fhaf is roughly proporfional fo fhe volume of B, up fo 
some small addifive error, called fhe discrepancy. 


Definition 20. Multi-dimensional discrepancy 

For integer s, put = [0,1)®. Given a sequence x = (xn)n=i, with Xn G F the discrepancy Dn{x) is 
defined as: 


F>n(x) = sup 
bgb 


^'^XB{xn)-yo\{B) \ , 


where XB{xn) is an indicator function which is 1 if Xn G B and 0 otherwise, and B is a non-empty 
family of Lebesgue-measurable subsets of I^. 


In this work, we shall define B as fhe sef of all s-producfs of infervals 


S 

[ui,Ui](modl) C [0,1)- 

i=l 


Offen in liferafure, one considers insfead fhe sfar-discrepancy D'f{x) which is defined by opfi- 
mizing over fhe sef of all infervals of fhe form [Oj 

The definifion of discrepancy nafurally admifs an interprefafion in ferms of probabilify: a 
sequence x = {a;„}„g[Ar] can be interprefed as a random variable x, fhaf fakes fhe value a;„, 
as n is fhe random variable n ~ U[N]. Saying fhaf Dm(x) < Dm means fhaf for any B G B 
fhe probabilify fhaf Xn is confained in B is equal fo vol(i3) -I- e, where |e| < Dm- In fhis work, 
deviafing somewhaf from sfandard ferminology, we shall offen refer fo fhe discrepancy of a 
distribution x on lengfh N, s-dimensional sequences. In fhis case Dm(x) will denofe 

Dm(x) = sup {\Px'-.x,n-^U[N](Xn G B) - VOl(B)|} , 
bgb 

''applying the union bound over all eigenvalues 
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Low-discrepancy sequences have much in common with random sampling, or the Monte- 
Carlo method, in the sense that they visit each cube a number of time that is roughly propor¬ 
tional to its volume, up to a small additive error. Yet, contrary to the Monte-Carlo method, such 
sequences are not random, but only appear to be random in the sense above. Arguably such 
sequences are most useful in fhe confexf of numerical infegrafion: insfead of compufing an in- 
fegral of a confinuous funcfion / over fhe s-dimensional cube [0,1)® one can replace if wifh fhe 
average value of fhe funcfion over fhe poinfs of a low-discrepancy sequence x = The 

well-known Koksma-Hlawka fheorem, fhen connecfs fhe approximation error as a funcfion of 
fhe discrepancy Dn{x) as follows: 


1 

N 


N 

- / 

i=i 


f{z)dz 


<V{f)-D*^{x), 


where V (/) is fhe bounded-variafion of / on [0,1)®. 

There are deferminisfic s-dimensional sequences x = {xi}f^i wifh discrepancy as low as 


Dn{x) <C- 


log^iV 

N 


and mafching lower-bounds (up fo consfanf factors) on fhe smallesf possible discrepancy are 
known for s = 1 [Nie92]. Hence, usually one considers low-discrepancy sequences fhaf are 
very long {N) compared fo fhe dimension (s). 

Af this point, it may be insightful fo consider an example of a low-discrepancy sequence: 
fhe well-known van-der Corpuf sequence [Nie92]: consider fhe binary expansion of a positive 
infeger 

b 

Vne[A = 2^] n = 

i^Q 

then the van-der Corput sequence x = is defined as: 


b 

VnG[A = 2'’] a;„ = ^a,2-*-^ 

i=0 


This sequence has 


D*n{x) < C 


\ogN 

N ■ 


We note fhaf fhe discrepancy upper-bound decays asympfofically like 0(1/A) (assuming 
small dimension s) whereas fhe Monfe-Carlo mefhod converges more slowly, behaving as 
0{1/\/N) - and hence fhese sequences are often preferred as a mefhod of numerical infegrafion 
fo Monfe-Carlo. They are also advanfageous compared wifh purely deferminisfic mefhods, like 
defining a fme-resolufion grid, because usually one can increase fhe lengfh of fhe quasi-random 
sequence, and reduce ifs discrepancy while making use of all previous poinfs of fhe sequence. 


4.2 Some basic facts 

We require a Lemma [2.5] due fo Niederreifer [Nie92]. 

Lemma 21. [Nie92] Small point-wise distance implies similar discrepancy 

Let xi,... ,X]S[, yi,..., j/AT denote two s-dimensional sequences for which \xn^i — Unf < O for all 
n G [A], i G [s]. Then the discrepancies of these sequences are related by: 

\Dn{xi, .. .,xn) - I?Ar(j/i,..., j/Ar)| < s • e. (1) 
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We prove an additional fact: 

Fact 22. Let x = be a distribution on sequences with discrepancy at most Dn{x), and let 

y = {yn}nelN] denote the fallowing sequence: 

yn=Xn + Zn, 

where z = {zn}n^[N] is some sequence chosen independently from x. Then 

DN{y) = Dn{x). 

Proof. For each S G B we have 

Py,n~U[N]{yn € S) = Px,z,n~U[N]{Xn~LZn GS)= Px,n'^U[N]{Xn G S zfa'P^n'^U[N]{Zn)dZn 

Jlo.iy 

By the discrepancy assumption on the sequence x we have that 

\/S G B \Px,n~U[N]iXn G S) - VOl(S')| < D]s{x). 

Therefore by fhe above 

Py,nr..U[N]{yn GS)= (vol(S') + e(z„)) • PzfazfadZn, {cfafal < Dn{x). 

J[0,1Y 

and hence by convexify: 

Py,n^U[N]{yn G S) = VOl(S') + £, |e| < Di^{x), 

which implies fhe claim by fhe definition of discrepancy 

□ 


4.3 The Good Seed Problem 

In our confexf we will be inferesfed in sequences x = {xn}n=i of fhe form 



where g G [fV]® is some s-dimensional vecfor, called fhe seed of fhe sequence. Specifically, fhe 
vector g /N will be fhe vecfor of eigenvalues of an n x n Hermifian mafrix A whose specfrum 
g = C{A) we would like to analyze. Since if is unreasonable fo assume fhaf fhe inpuf mafrix has 
a specfrum fhaf is a good seed, we would like fo find a perfurbafion of fhe mafrix A' = A-\- S 
such fhaf g' — C{A') has a corresponding sequence, defined as above, wifh low-discrepancy 
Niederreifer has shown [Nie92] fhaf if g is sampled uniformly on [TV]'^ fhen if is a good seed 
wifh high probabilify: 

Lemma 23. Let s,N bean integers and g ~ [/([iV]®). Then 


P 


(^Dn{x) < 


log‘‘N'\ 

) 


> 1 - 1/N. 


For our applicafion we require fhaf N = poly(n), and s = 2, in which case fhe above 
discrepancy is sufficienfly low for our purposes. Yef, since if requires fhe normalized seed g/N 
fo be essenfially uniform on [0,1)", if implies fhaf fhe corresponding mafrix perfurbafion S 
added fo A musf be very sfrong - fhereby loosing all cormecfion fo fhe inpuf mafrix. 


14 



4.4 Finding Reasonably-Good Seeds Locally 

To bridge the gap between weak-perturbation and low-discrepancy we show a new lemma, 
which may be of independent interest: it allows to trade-off fhe exfenf fo which g is random, 
and fhe discrepancy of fhe sequence generafed by g. Specifically, we will show fhaf if g/N is 
uniform on cubes of much smaller side-lengfh, i.e. af leasf 1 / VN, fhen fhe resulfing sequence 
has discrepancy 0{log^N/'/N). This is fhe subjecf of fhe following lemma: 

Lemma 24. We are given integer N, with prime divisor M = Q{N°‘)for some constant a > 0, and an 
integer s. Let g = (gi,... ,gs) G N^, such that each coordinate gt is independently chosen uniformly 
on some interval It C [N] of size M. Let x = x{g) = {xn}n=i be the following s-dimensional sequence 
of length N corresponding to residuals of g: 

N /■ 


Then 


Pg (^Dn{x) < 21og®(M)/yM) > 1 - 1/Vm. 


Proof. For infegers P, s puf C* (P) as fhe sef of all vectors in Z® wifh entries in [—P/2, P/2) n Z, 
excluding the all-zero vector. Following Niederreiter [Nie92] we define for each h GiP 


rih) = Wr{hi), ( 2 ) 

2=1 


where r{hi) = max(l, \hi\). For g — (^i,..., we denote: 


Ri9,P)= E (3) 

h-g—0{modP) ,h^C* (P) 


By Theorem [5.10] of [Nie92] when each gi is randomly chosen on fhe entire interval \P], for 
prime P, fhen 

E,[i?(g,P)]=o(^i5ffl). (4) 

Since R{g, P) > 0 for all vectors g, fhen 

Pg{R{g,P)>^^^^^<l/y/P. (5) 

or 

Pg [nig, P) < > 1 - i/v^- (6) 

Lef us use fhe above equation fo upper-bound fhe discrepancy of S{g). Recall fhaf M is a prime 
divisor of N, wifh M = <d{N°‘). We firsf observe fhaf: 

R{g,N)= E rih)-^< (7) 

h-g=0{modN),hGC*(N) 


E E ^ 

h-g=0{modM),heC*{M) h-g=0{niodM),hGC^(N),3i, S.t \hi\>M 

We note fhat fhe second ferm is at most 


s 

M 


E ^ 

;iGZ=-i 


s-log®-i(iV) 

M 


( 8 ) 


(9) 
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Regarding the first term, h ■ g = O(modM) if and only if h ■ (g{m.odM)) = O(modM), and since 
each gi is uniform on some interval of size M, fhen gi{modM) is uniform on \M]. Since M is 
prime, we can apply equation (6) to the first term with P — M. Using this equation and the 
upper-bound on the second term we have: 

Pg N) < > 1 - 1/V^. (10) 

According to Theorem 5.6 of [Nie92], the discrepancy of fhe sequence {gn/N}^~Q is upper- 
bounded by: 

D^{S{g))<j^+2-^-R{g,N). ( 11 ) 

Plugging equation (10) into equation (11) implies: 


(^D^{S{g)) < > 1 - I/Tm. 


( 12 ) 


□ 


5 A Filtering Algorithm 

In this section we provide the specification of fhe filfering algorifhm, which is fhe main compu- 
fafional black box of our algorifhm. This algorithm accepts an integer m that separates the z-th 
eigenvalue of a Hermifian mafrix A and compufes an approximafion for fhe z-fh eigenvector, 
with high probability: 

Algorithm 25. Filter(A, m, d) 

1. Compute parameters: 

p = 247z^|’ln(l/(5)],^ = S‘^j{2pm). 

2. Sample random unit vector: 

Sample a standard complex Gaussian vector v, set wq = z;/||z;||. 

3. Approximate matrix exponent: 

Compute a C Taylor-series approximation of e'’^, denoted by U. 

4. Raise to power: 

Compute U'^ by repeated squaring. 

5. Generate matrix polynomial: 

Compute B = ^ ^ by repeated squaring. 

6. Filter: 

Compute zc = 

7. Decide: 

Set z = A ■ w, io = argmaxig[„] \wi \ and compute c = Zigfwi„. If 

IIA • zo — c • zoll < SJyhz 

return w, and otherwise reject. 
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We now show that if the algorithm is provided with an integer m that separates the fc-th eigen¬ 
value of A in fhe sense defined in Definition 6 , fhen fhe outpuf is close fo fhe fc-fh eigenvecfor 
of A. 


Theorem 26. Let S < and a = ^ln(l/^). We are given annxn Hermitian matrix A with 
eigenvalues and a corresponding orthonormal eigenbasis {vi}i^yn]- Additionally, we are 

provided an integer m that separates k in A, w.r.t. Bout, In the sense of Definition 5. Let 

w = Filter{A,m,S). Then 

P (llw — r;fe|| < (5) > 1 — 3n“^. 

The algorithm has boolean complexity at most 2cn‘^ ■ log(2p^m^/5^), and runs in parallel time 
0 (log^(n)). 


Proof. Lef {rf}£g[„] denofe fhe sef of eigenvalues of U. Since C/ is a pol 5 momial in A (fruncafed 
Taylor series) fhen {ve}((=[n] is also an orfhonormal basis for [/. Since in addition |l^-e*^|| < C 
fhen 

V£e[n] (13) 


Lef w' = B ■ wq and denofe 


Wo = ^ PtVi, and w' = ^ acvc. 

ieln] £e[n] 

Since A, U share fhe same basis of eigenvecfors, fhen by fhe detinifion of fhe mafrix B fhe 
coefticienfs /3i are relafed by: 




l + T^ 


2p 


2 


So by Equation 13 




2 ^ 


I _|_ gimXi 


2p 


— 2pm(!^. 


Since m separafes k fhen {mAfc} G Bin, and for all £ ^ fc we have {toA^} f. Bout- Thus, for i = k: 






2 ^ 


-I- cos(l/2n-\/ln(l/^)) 


— 2pmC, 


Using Claim 30 


> 1 - 


1 


4n^ ln(l/J) 

On fhe ofher hand, for all f 7 ^ /c we have: 

1 2 2 -(- girnXe 


-2pmC>—,. 


(14) 


m 


< 


-|- 2pm(^. 


SO since m separafes k fhen: 


<1(1-1- cos(l/2n))|^ -I- 2pm( 
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(15) 


Using Claim 30 

< (1 - l/12n2)24""in(i/'5) 2pmC, < e-2i"(i/-5) + 2pmC, < 25'^. 


By Fact 29 for any e = 1 /poly(n) there exists a constant c > 0 such that 

P(Vh j \l3j\ < c\|3i\^J\n{l/£)/£) > 1 - ine. 
Choose £ = n~'^. Then by Equations 14 and 15: 


p(y£^k 




< c‘‘{2S^) ■ (2e°) • 41nn • > 1 — 3n 


and so for S < n there exists p G C, \ri\ = 1 such that 

2 

P -Vk < 


1 


m\ 


< 16c^n® Inn^^e® < S. 


i/fc 


for sufficiently large n. Using Claim 28 we conclude that w.p. at least 1 — 3n~^ over choices 
wq, the criterion is met and the algorithm returns a vector w = i(;'/||iy'|| satisfying the equation 
above. 

Arithmetic run-time: The approximation of by U requires, using Fact 31 a time at most 

cn“log(l/(^) = cn‘^ ■ \og(2pm/S^). 

Next, the repeated powering of C/ to a power m requires time at most: cn‘^[log(TO)] and the 
repeated powering of B to the power p requires time at most: cn'^ [log (p)] Hence the total 
complexity is : 4cn‘^ • \og{pm/6^). 

Depth complexity: Each matrix product can be carried out in depth log(n). Each of steps 3 to 6 
involves at most log(TO) + log(p) sequential matrix multiplications. The sorting algorithm can 
be computed in depth 0(log(n)) by Fact 11. Hence the depth complexity of the entire circuit is 
at most log(n) • (log(m) + log(p)) + 0(log(n)) = 0(log^(n)). 


We conclude the proof of the theorem by showing stability: 

Claim 27. Under the assumption of Theorem 26 the algorithm is log-stable. 

Proof: Consider the arithmetic operations involved in computing the filtering algorithm: 

1. Generating an approximation [/ of as a truncated Taylor series. 

2. Raising U to a power m G [M]. 

3. Computing ((/ +U)/2)^’. 

4. Normalizing ilwo/||t?rcoll- 

Consider an arithmetic circuit C implementing the above, and the circuit D = D{C,t) - 
the discretization of C to t bits of precision modeled as follows: after each arithmetic step, 
the result is rounded to the nearest value of 2“*. Consider all steps except division. A is 5- 
separated so in particular ||A|| < 1. Thus, whenever we multiply two matrices at any of the 
steps above both have norm at most 1. Hence, at each rounding step the error is increased by 
at most yTi2“h Finally, considering the final division step, we observe that since m separates 
k, then by Equation 14 we have ||Bwo|| > l/(2e®). This implies that the total error is at most 
y/n{p + M) ■ 2“* • 2e®. Since M,p are both pol 5 momial in n then for any S = l/poly(n) the error 
is at most S for some t = 0(log(l/S)). 

□ 
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5.1 Supporting Claims 

Claim 28. Suppose that ||w — ti|| < 5 for some unit eigenvector v of A, and S < 1/4. Let z = A - w, 
and iQ denote io = argmaxig[„] |wi|. Let c = Zi„lwi„. Then 


Proof. Write w 
Therefore 

So 


\\A ■ w — c ■ w\\ < 3S^/n. 

V + £, ||£|| < 6. Since |jA|| < 1 then z = Aw = Xy ■ v + £', where ||£'|| < S. 

ZiQ = At, • Vig + £ig = At, • {wig — £ig ) T £^g , 

Zig = XyWig + £”, \£''\ < 2S. 


Since IwigI > l/-ynthen: 

Hence 


c — —^ — A„ + (/, l/^l < 26y/n. 

W,g 


||t4 • w - c • w|| < \\Xy ■ V + £' - (At, + C) • (v + f )|| = \\£' - (Xy + ()£ - Cull 

which by the triangle inequality is at most 

5 + (1 + 2d^/n)6 + 26y/n < 3S^/n, 


for sufficiently large n. 


□ 


Fact 29. Random unit vectors have well-balanced entries 

Let be some orthonormal basis o/C”, 0 < e = l/poly(n), and u G C” a uniformly random 

complex unit vector. For any i G [n] let Ui = |(u, Ui)|. For any e = 1/poly(n) there exists a number 
Cl > 0 independent of n, such that 

P (yi,j Inil/luil < ci^J\ii{l/£)/£^ > 1 - 3n£. 

Proof. Sample a random unit vector by sampling a standard Gaussian vector and normalizing 
to unity Let x = (xi,..., x„) where Xi ~ p(0,1), and all are i.i.d. Since the Gaussian measure 
is invariant under conjugation by a unitary matrix we can assume w.l.o.g. that Vi = Ci for each 
i, in which case the afs are i.i.d. at ~ p(0,1). By the Gaussian measure for every C 2 > 0 there 
exists Cl > 0 such that 

P(|ai| > ciV• Inn) < n~‘^^. 

Therefore, by the union bound 

P(Vt |ai| > ciV• Inn) < 

On the other hand, for any ai, 0 < e < 1/n, we have 

P(|ai| < e) < 2e. 


So by the union bound 

P(Vi G [n] \ai\ < e) < 2n£. 

Set C 2 = ln„(l/£). Since £ = n~^ for some constant k, then C 2 = fc is also a constant. Then there 
exists a constant ci such that: 

P j \ai\/\aj\> ciV\nm/£^ < 2n£ + = 3n£, 

Since normalization does not change these ratios, then the maximal ratio between the compo¬ 
nents of X, is the same as the maximal ratio between the components of x/||x||. This implies 
the proof. □ 
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Claim 30. V6» G [-0.01,0.01] 1 - ^ < < 1 - f • 


Proof. Follows by truncating the Taylor series of cos(a;) to second order. 


□ 


Fact 31. Efficient approximation of exponentiated matrix 

Given a Hermitian n x n matrix A, ||^|| < 1, and error parameter e > 0, a Taylor approximation of 
denoted by Ua can be computed in time cn‘^log(l/£) and satisfies 


- Ua 


< e. 


Proof. Put s = log(l/^) and consider the Taylor approximation of up to the first s terms: 

s-l 


Ua:=J2 


(tAl 

7,1 


m—0 


Then the approximation error can be bounded as: 


- Ua 


I ^ 'YY) I ! ^ YU ^ 


OO - 

< V A < 2"® < 

m l 


(16) 


where we have used the fact that ||£l|| < 1. The complexity of fhe approximation is comprised 
of log(l/£) mafrix producfs for a fofal of cn‘^log(l/£). □ 


6 Sampling Separating Integers 

In fhis section we show our main fechnical fool: which is fhaf perfurbing a (5-separafed Hermi- 
fian mafrix Ahy a Gaussian mafrix of a carefully calibrafed variance, resulfs in a corresponding 
sequence of residuals S'(gI) having low-discrepancy, af leasf for pair-wise variables, which in 
fum implies fhaf we can separafe each eigenvalue of A almosf uniformly: 


Theorem 32. Let A be a 6-separated n x n PSD matrix, E GUE, ( < and 

4 < a < n. For any M > we have: 

Wk G [n] Pe,m~u[M] ( cn separates kin A + ■ £ w.r.t. Bin{a),Bout ) > l/(5an) 


6.1 Additive Perturbation 

Definition 33. {a, s)-normal vector 

Letv= {vi,... ,Vn) be a vector of n random variables, vissaidto be {a,e)-normal,ifeachvi ^Xi + yi 
with Xi ~ p(Ai, a'^)for some G K, xi,..., a;„ are independent, and lyf < as, for all i. 

By our definitions above, Gaussian perturbation of a mafrix wifh well-separafed eigenval¬ 
ues resulfs in a (a, £)-normal vecfor as follows: 

Fact 34. Perturbation of well-separated matrices 

Let A be an n X n aminSeparated Hermitian matrix with eigenvalues Ai > A 2 ... > A„, where 
oimin > £/ cind £ > 0 is some constant. Let £ be GUE, and A' = A + s^ ■ £, where L>2. Then w.p. at 
least 1 — n- 2“” the vector of eigenvalues of A' {X[,..., A(,) is cne^~^)-normal, for some constant 
c > 0. 
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Proof. We invoke Corollary 16 choosing e as and 5 as e. We get: 

Vi G [n] A-= (1 — a) •/i(Ai, + a • I? + Ci, |Ci| < 16cn • 0 < a < 2“". 

Choosing a = implies the each A' w.p. at least 1 — 2“" is sampled according to: 

A' = + Ci, ICil < • 16cn£^“h 

Taking the union bound over all i G [n] then implies the proof. □ 


6.2 Moderately Low-Discrepancy 


Lemma 35. Low-discrepancy sequence from almost normal vectors 

Let B > t), and v = {vi,..., Vn) be some (ct, e)-normal vector, for a = n~^ ,e < There exists 

M < such that for any S = {ii,... ,is} C [n] , j^l = s the distribution on s-dimensional 

sequence of length M: 

W = {{{m -Vi,} 


satisfies 


Dm{Vs) < 41og®(n) • n 


Proof. Let P be the minimal prime at least , and put M = P®. By Bertrand's postulate, for 
sufficienfly large n we have fhaf M = P^ < For any z G [0,1) lef be fhe 

number closes! fo z in fhe grid m/M, m G [M]. 


Removal of non-independent component. Since v is (cr, e)-normal fhen Vi = + Yi, 

where Xi ~ \Yi\ < etr, and fhe Xfs are independenf. Lef denofe fhe sequence 

generafed by faking only the X component of fhe seed vecfor v, i.e.: 

ys = (17) 

Fact 36. 

Dm{Vs) < DMiV^) + s • 

Proof. Consider the r.v.'s Xi,Yi. By our assumption 

Vi G [n] \Yi\ < ae = . (18) 


Thus: 

Vm G [M], i G [n] \{mvi} — {mXi}\ < m ■ (19) 

By Lemma 21, we can conclude that the discrepancy of our fargef sequence Vs follows tighfly 
fhe discrepancy of : 

DM{Vs)<DM{Vi^) + s-n-°-^^ ( 20 ) 

□ 


Reducing Gaussian measure to uniform measure Consider fhe vecfor derived by frun- 
cafing each coordinafe of fhe vecfor {Xi^^ ,■■■, Xi^) fo fhe nearesf point on the M-grid: 

X^ = (Xf,... ) = (IMX,J/M ,..., IMX,J/M). 

Consider the discrepancy of fhe disfribufion on s-dimensional sequences formed by faking 
infeger mulfiples of X^. We claim: 
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Fact 37. 


P,; {pMiVg’^) < log'’(n) • n > 1 - 3n 

Proof. In Fact 40 choose as parameter m = We get that w.p. at least 1 — 2n^/m = 

1 — each Xi samples a convex mixture of variables {wj}j^[m] where 

Wj ~ U{Ij), \Ij \ = aim = ( 21 ) 

Hence, w.p. at least 1 — 2n“°-^^ for all i G [n], fhe variable M ■ {Xf^} is a convex mixfure of 
uniform random variables on intervals Ij C [M], where 

|/.I > ^ > „i-5n . „-i.2B-2 > (22) 

m 

We apply Lemma 24 to the sequence of residuals of infeger mulfiples, wifh fhe seed X^^: 

^ ({mWf } ,..., {mXf . (23) 

The lemma requires fhaf each variable be disfribufed as: 

MXf^ -- t/[X], 

where I is some inferval of \M], for infeger M > 1 safisfying: 

\I\ > P, P prime , P > N°‘, a > 0. 

By our choice of paramefers M has a prime divisor P equal fo = P. Hence, by Equafion 
22 we can safisfy fhe assumpfion of fhe lemma by choosing fhe paramefers N, M, a for Lemma 
24 as follows: 


N = M,M = P,a = 0.2. (24) 

Hence, by Lemma 24, and accounfing for fhe Gaussian-fo-uniform approximafion error we gef: 

Pv < 21og®(n) • > 1 - > 1 - . (25) 

□ 


Treating the residual w.r.t. the M-grid Define: fhe fruncafion error 

V^e[s] r, :=X,-Xf. 

In Facf 36 we analyzed fhe error Yi whose magnifude is negligible even w.r.f. 1/M, and can 
fhus be disregarded for any elemenf of fhe sequence Vs- Unlike fhis, fhe residual error r; cannof 
be disregarded because when mulfiplied by infegers uniformly in [M] if assumes magnifude 
U(l). Thus, if requires a differenf freafmenf. 

Corollary 38. 

P. {Dm{V^) < 21og*(n) • > 1 - 

Proof. Express fhe z-fh elemenf of fhe sequence using r^: 

Vz G [s] {X.,-m} = {{Xf^ + ri)-m]={{mXf^} + {mri}} (26) 
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Let E denote the event in which Xi is sampled according to Wj ^ U\Ij\ where wj is at distance 
at least 1/M from either one of fhe edges ofij. Condifioned on E, fhe random variables and 
X:^ are independenf for all i G [s] so condifioned on E we also have 

Vi e [s] {mX^}, {mri} are independenf. 

is a disfribufion on sequences formed by sampling fhe inifial seed {Xf^}, and is a 
disfribufion on sequences formed by adding fo a random sequence formed by sampling 

fhe independenfly seed and fhen generafing fhe lengfh-M sequence 

{(mri,... ,TOrs)}mG[M]- 

Then since condifioned on E {mri} are independenf for all i, fhen by Facf 22 we have: 

DMiVi^\E)=DMiV^^^) 


and so by Facf 37 

P„ {DMiVi^lE) < log*(n) • > 1 - 

By Equafion 22 fhe probabilify of E is af leasf: 

PviE) > 1 - |I,|/(2M) > 1 - M°-V(2M) > 1 - n"'®. 


Thus: 


Pv {Dm{V^) < log"(n) • > 1 - - P(£;) > 1 - 


(27) 


(28) 

□ 


Conclusion of proof: 

By Corollary 38 we have 

P. {Dm{V^) < 21og*(n) • > 1 - 

and by Facf 36 we have 

Dm{Vs) < Dm{V^) + s • 


Thus by fhe union bound: 

Pv {Dm{Vs) < 21og"(n) • + s ■ > 1 - dn’ 


fhus: 

P^ {Dm{Vs) < 31og®(n) • > 1 - 4n“°-^^ 

Hence for all buf a measure 4n“° of sampled vecfors v, fhe resulfing sequence has discrep¬ 
ancy af mosf 31og®(n)n“°'^^. Since fhe discrepancy measures fhe worsf-case addifive error for 
any sef fhis implies fhaf: 

Dm{Vs) < 31og®(n)n-° i^ -h4n-° i^ < 41og"(n)n-°-i^ 


□ 
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6.3 Approximate Pairwise Independence 

Lemma 39. Let A = (Ai,..., A„) G [0,1)” and M a positive integer that satisfy: 
yii^j Dm {{{mXi,mXj)}me[M]) < C, 

Let 4 < a < n. For each k G [n] w.p. at least l/(5an) over choices ofm ^ U[M] the sampled sequence 
m separates k w.r.t. Bifya), Bout- 

Proof. Let Ei denote the following event: 

Ei .— Xi G Bj^ji A yj f i Xj Bout 

We want to show that 

Vf G [n] PiEi) > 

oan 

Let s denote the number of x/s not in Bout- 

~ |{7 3 i ^ 

Then under fhis nofafion we have: 


P{Ei) = P(s = n - 1 A Xi G Bin)- 

Consider fhe conditional expecfafion: 

E [s|Xi G Biu\ 

By linearify of expecfafion: 

E [s|Xi G -Sin] ~ 'y ] P [^j f- Dout\^i ^ Sin] ■ 
f A* 


(29) 


Considering each summand separafely: 

^ Sonijxi G Sin) — 


^ Dout X Xi G Sin) 


P(Xi G Sin) 

Using fhe pairwise discrepancy assumption, fhe above is af mosf: 

(l-|Sn„i|).|S„l+C ,0 1,0. 1,0./ 0-51 

-——- - - < 1 - |Sont| + 2Qan = — + 2q;C’^ < - 

|Sjn|-C 2n n 

and so by Equafion 30 

E [sjXi G Sin] = (n - 1) • P [Xj f Bout\x^ G S^n] < 0.51. 

The variable sjxi G Sin accepfs only infegral values, and by Markov's inequalify: 

P(s > IjXi G Sin) < 0.51 

Therefore 

P(s = Ojxj G Sin) > 0.49. 

Using again fhe 1-dimensional discrepancy we have 

P(xi G Sin) > -C > y—• 

an 2an 

Subsfifufing fhe lasf fwo inequalities into Equation 29 yields: 

P(Si) > 0.49 • ^ 

Zan ban 


(30) 


□ 
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6.4 Proof of Theorem 32 

Proof. By assumption A is ^-separated and < min{n“®°, Consider the perturbed matrix 

A' = A + C£. 

Choose L = 13 and e = Then by Fact 34 w.p. at least 1 — n2“" all eigenvalues of A' are 
(cr, e)-normal with parameters 


a = C < n"™,e < < C°'® = ® 

where the last inequality follows because ( < We assume fhaf fhis is fhe case and accounf 
for fhe negligible error af fhe end. Sef now B = log„(l/cr). Then fhe eigenvalues of A' are (cr, e)- 
normal for cr = n~^ an e < . Since in addifion 4 < a < n fhen by Lemma 35 fhere exisfs 

aninfeger M < safisfying: 

VS' C [n] DM{{mXs}) < 41og'*(n)n“^ < n~^, (31) 

for sufficienfly large n. Hence, by Lemma 39 a random m ~ U[M] separafes fhe k-ih eigenvalue 
of ^ + S w.r.f. Bin{a), Bout w.p af leasf l/(5an). 

□ 


6.5 Technical approximations 

Fact 40. Approximating a Gaussian by a convex sum of uniform distributions 

Let g = {gi,... ,gn) be standard Gaussian vector g G K". Then g is equal to a convex combination 
oftivo distributions as follows: (1 — p)'Du + p ■ Vy, where Vu is the n-fold distribution of 

independent variables zi,... ,Zn, where each Zi, \zi\ < y/n, and Zi is the convex sum Zi = tjWj 
ofm > 2rf Lid. uniformly distributed variables, with Wj ~ U[Ij], with \Ij\ = Ijm, andp < 2vf /m. 

Proof. Partition the interval [—vTi/2, yTi/2] into m ■ ^/n equal intervals Ij, each of size 1/m, 
and lef pj denofe fhe poinf of minimal absolufe value in fhe j-fh interval. Set tj = ^{pj). 
Consider the real Gaussian PDF $(a;) = For any pair of neighboring infervals 

Pk,Pk+i, \Pk+i \ > |pfc|, wehave 


\<t>{pj) - $(^^+1)1 < • 



e-(Pj+i/"i)V2 

e-tPifn 


< 1 _ e-Pr/"* < 1 - < vS^/m. 


(32) 


Then by fhe above: 


Vj G [m] max{|tj — \tj — tj+i\} < y/njm, 

Lef Zi : K I—[0,1] denofe fhe following funcfion, which is a sum of uniform disfribufions on 
fhe infervals Ij: 

rriy/n 

Zi = tjU[Ij]. 
i=i 

Consider fhe /i-disfance between gi and zp 



Zi(a;)| dx 



dx 
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By Equation 32, and by bounding the Gaussian tail at we establish the upper bound: 


^ m^/n 

< — ^ max{|tj — \tj — tj-iW + 2 • e”" < ^njm. 

7Tl 

On the other hand, by definition of fhe poinfs tj, we have 

gi{x) — Zi{x) >0, Vx S K. (33) 

Hence each gi may be wriffen as a convex combinafion 


g^ = (l-p) ■ Zi+p-y^, z^{x) 


Ztjx) 

J^Zi{x)dx 


P < 2n/m. 


Since fhe variables gi are independent then the n-fold distribution of ,..., gn, can be written 
as a convex combination of fhe of fhe n-fh fold disfribufion of such i.i.d. variables zi,..., Zn, 
and anofher disfribufion, occurring, by fhe union bound, w.p. af mosf n ■ p < 2n^/m. □ 


7 Parallel Algorithm for ASD 

The algorifhm Filter(,4, m, i5) described in Section 5 is given an infeger m fhaf separafes fhe 
f-fh eigenvalue, and refurns an approximation for fhe i-fh eigenvecfor. In fhis section, we use 
fhis algorithm in a black-box fashion and design a Las-Vegas algorifhm for computing fhe full 
ASD of a mafrix. Essentially, if amounfs fo rurming sufficiently many copies of fhe tittering 
algorithm in parallel so that all eigenvectors are (coupon) collected with high probability. 


Algorithm 41. Input: n x n Hermitian matrix A, parameter 5 <l/n. 

1. Compute parameters: 

B = min{5, i?*(5/(3-\/n))},i5^ = (min{i5, i?})^^/4 

a = = (max{B“^^,n“®°})^®,T = 60nQ;log(n) 

= [1 - l/(2n),-1 -I- l/(2n)] and = [-l/(na), l/(na)], 

2. Perturb: Ai = A + ^{8'Y + <5^/(9n) • E, where E is CUE. 

3. Run T parallel processes of the following procedure 

(a) Sample m ~ U[M] 

(b) Run Filter (Ai,m, S'). 

4. For vector w = Wk sampled at process i S [T], compute z = A ■ w, = argmax^gj^j \wi\ 
and Afc = Zi^jw^^. 

5. Sort the values Xp assume Ai < ... < A 7 -. Initialize: 7 = Ai, I? = $. Iterate over all 
i = 1,... ,T. At each step i: if \^ — Xi\ > B/A then add V VU {wi}, and set 7 = A^. 


We now state our main theorem: 
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Theorem 42. For any Hermitian matrix 0 ^ ^ ^ 0.9/, and 5 < n there exists an 
algorithm computing ASD(A, (5), in boolean complexity 0(n‘^^^). The algorithm is log-stable. 


Proof. 

Correctness: 

Let £i , £2 be independent GUE matrices, and set 

Aq = A + (S/(3\/n)) ■ £ 1 . 

By Lemma 18 and definition 19 the parameter B < B* above satisfies: 

P (Aq is B separafed) > P (Aq is B* separafed) > 0.99. 

Assume fhaf fhis is fhe case, and accounf for fhe error af fhe end. Consider now a small-scale 
perturbafion of Aq: 


Ai — Aq S' • £ 2 . 

Since £i ,£2 are independent the matrix Ai is a perturbation of A as follows: 


(34) 


Ai = A + ^((502 + (<5/(3v^))2.£, (35) 

for GUE mafrix £. Hence fhe mafrix Ai used by fhe algorifhm sfarfing af sfep 3 is a 5'- 
perturbafion of a il-separafed mafrix Aq. 

We now invoke Theorem 32 w.r.f. Aq. Aq is separafed wifh paramefer B, and perfurbed by 
GUE wifh sfandard deviafion ( = S' < B^^. We also have S < which means in particular 
fhaf fhe paramefer ( used in fhe sfafemenf of Theorem 32 safisfies ( < n-^°. Einally we have 
2 < a = 0(-y/ln(u)) = o{n). Therefore: 

Vk G [n] Pe2,m~u[M] {m separafes A: in Ai = Aq -I- S' £2 w.r.t. Bin{a), Bout) > l/(5an) 

Lef denofe an orfhonormal basis for Ai, and lef Xi denofe fhe eigenvalues of Ai. 

Condifioned on sampling Ai, m such fhaf m separafes k in Ai, we invoke Theorem 26 for 
w = Filter(A, m, S'). By our choice of paramefers we have fhaf S' < and a = 

Hence, fhe condifions of fhe fheorem are mef and so fhe oufpuf vector w safisfies: 

Pu,o (Ik - < iJ') > 1 - 

Thus, by fhe union bound we have 

Vfc G [n] Pe2,m,wo (Ik - I'fcll < 3') > G-- 

ona bna 

Therefore, by fhe coupon collecfor's bound fhe probabilify fhaf T = 60nQ;log(n) parallel copies 
sample an approximation for each k G [n] is at least 

1 - 6/60-0.01 = 0.89 

Assume that this is the case. Let denote the output set of vectors from all T parallel 

copies. We now wanf fo show fhaf each oufpuf vector is a valid approximate eigen-vecfor. 
Each oufpuf vecfor safisfies by fhe stopping criferion: 


3c 


Xi ■ Wi- Ai ■ Wi 


< 3S'y/n < R/200, 
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where the last inequality is for sufficiently large n, and Xi is the algorithm's approximation of 
fhe i-fh eigenvalue. By Facf 15 since is S-separated fhen 

P(^i is B/2 separafed) > 1 — 2“". (36) 

Assuming fhis is fhe case, 5'y/ri is negligibly smaller fhan fhe separafion of Ai - so by Facf 44 
fhe value Xi above musf safisfy: 


3k G [n] Xi — Xk < B/10. 
and fhere exisfs some eigenvector Vk of Xk such fhaf 

\\wi - rifcll < B/30. 


This implies by fhe triangle inequality that for any j ^ i, I ^ k for which 
have 


A,- - Xi 


< B/10 we 


Xi A, 


> B/2-2B/10> B/A. 


hence Xj and Xj are classified fo differenf eigenvalue bins af sfep 5. We conclude fhaf fhe sef of 
values {Xi}]Li safisfies fhaf 


• For every k G [n] fhere exisfs j gT such fhaf: 


Xk A, 


<5/10 and there exists a sampled 


vector Wj, and a unit eigenvector Vk of Xk such fhaf ||wj — ffc|| < 5/30. 


• For every j G T fhere exisfs a unique k G [n] such fhaf Xk — Xj 
Therefore, af fhe end of sfep 5 we have |I?| = n and 

Wk G [n] BwGV ||w — ffell < 5/10. 


< 5/10 


Hence, wifh probability at least 0.89 — 2 " > 0.85 the algorithm returns a database V = 
{rci,..., Wn}, such that there exists an orthonormal basis {z;i,..., Vn} of Ai for which 


Vfc G [n] — Wk\\ < 5/10 < (5/10. 

In fhaf case V is an ASD(Ai, (5/10). On fhe ofher and, by Equation 35 we have that w.p. 1 — 2“" 
Ai satisfies: 

P-Aif <((5/3)^ 

Therefore w.p. af leasf 0.84 fhe dafabase TX is also ASD(A, y/ (5/10)^ + (5/3)^) which is in par- 
ficular an ASD(A, S). 

Run time: 

By Theorem 26 the parallel time of each copy of Filter(A, m, i5') is af mosf 0(log^(n)), and 
arifhmefic complexify 0{n‘^). Hence, fhis is fhe depfh complexify of all T parallel copies of 
fhe circuif for Filter(A, m). The sorting step 5 can be implemented in depth 0(log(n)) by Fact 
11, using 0(n) processors. Hence the complete algorithm runs in parallel time 0(log^(n)) in 
arithmetic complexity 0(n“+^). 

Finally, we show that the entire algorithm can be implemented in log(n)-precision: 

Proposition 43. The algorithm 41 is log-stable. In particular it runs in boolean complexity 0(n“+^). 

Proof. Each parallel process is log-stable by Claim 27. Hence, to show that the entire al¬ 
gorithm is log-stable, it is sufficient to show that for a Gaussian perfurbafion fruncafed fo 
t = 0(log(n)) bifs of precision, fhe sfafemenf of Theorem 32 sfill applies - namely a uniformly 
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random m ^ U[M] separates each eigenvalue, with high probability, for the parameter M set 
in the algorithm. To do that, write 

£t = £ + D, 

where £ is the non-truncated GUE matrix, and D is some matrix sampled from a disfribufion 
such fhaf \Dij\ < for all i,j. Following Equation 34 define 

^ = d.g -\- £ -\- D = Ai + D, 

which is fhen used by fhe parallel process of Filter(A^*\ to, S'). By Equation 36 w.p. af leasf 
1 — fhe mafrix Ai is separafed wifh paramefer af leasf B/2. Hence, for sufficienfly large 
t = 0(log(n)), such fhaf 2* > B and using again Facf 15 we have for each i G \n\. 




< 2y/n-2~ 


Hence by Lemma 21 if 

{{mXi^{Ai)}, {toA*2 {Ai)})^^[m], 

is a 2 -dimensional sequence wifh discrepancy Dm, then the discrepancy of fhe noisy sequence: 
({toA.,(4 ‘))},...,{toA,,(4 ‘))})„,, [M] is af mosf Dm + 2 ■ M ■ ■ 2 *. This implies fhaf fhe 

pairwise discrepancy error increases by an addifive error af mosf 2MyTi2“‘ af Equafion 31 in 
fhe proof of fheorem 32. Hence, fhere exisfs a choice of f = 0(log(n)) so fhaf fhe discrepancy 
of fhe above sequence is af mosf n“^, fhereby satisfying Equafion 31. Hence, fhe sfafemenf of 
Theorem 32 holds for £ = £t, for t = 0(log(n)). 

This implies fhaf fhe algorifhm can be compufed in arifhmefic operation, where 

each operafion is a binary operation befween fwo regisfers of size 0(log(n)). Hence, fhe 
boolean complexify of fhe algorifhm is 0(n‘^+^) • log(n) = □ 

□ 


7.1 Supporting Claims 

Fact 44. Let A be a S-separated matrix, w a unit vector satisfying 

llAw-cwll <B,B< ^VlOO, 
for some |c| < 1. Then there exists A e C{A) such that 

|c- A| < (5/10, 

and there exists some unit eigenvector v of \ such that 

Ik ~ '*^11 ^ (5/30. 

Proof. Without loss of generalify we assume w = aivi + a 2 V 2 where vi,V 2 are fwo disfincf 
elemenfs of some orfhonormal basis of A, corresponding fo eigenvalues Ai, A 2 . Then 

II Aw - cw|k = Ikiai(c - Ai) -I- V 2 a 2 {c - A 2 )|k < 

By fhe facf fhaf vi , V 2 are orfhogonal: 

|ai4|c-Aip + |a2nc-A2p<B2. (37) 

Since ||w|| = 1 fhen \aif |Q! 2 p = 1, so for af leasf one i G {1,2} we have 

|c- A4 < 2 Bk 
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Suppose w.l.o.g. that i = 1. Since A is 5 separated then |Ai — A 2 I >5. Then together with the 
triangle inequality we have 

IA 2 - c| > IA 2 - All - |c- All >6-V2B > 0.95. 

Hence by Equation 37 above: 


|a2p < BV(0.95)2 < 5 V 1 OOO 

This implies that for some unit eigenvector ?; of Ai we have 

||w — w|| < 5/30, 

Using fhis inequalify back in fhe assumption, and since || A|| < 1, |c| < 1 fhen 

\\cw — Aw\\ = 110 ( 1 ; + S) — A - (v + S)\\ < B, |j£|| < 5/30 


and since Av = Xiv fhen 

||cw-Ai?;|| < B + 25/30, 

and since ||z;|| = 1 we gef: 

|c-Ai| < B + 25/30 < 5/10. 

□ 


8 Numerical Experiments 

In fhe work above we have ouflined a new fheorefical approach for solving fhe eigen-problem 
of Hermifian mafrices. Given fhe immense practical imporfance of fhis problem, we now pro¬ 
vide in addition some numerical evidence abouf fhe performance of fhe algorifhm. We run 
fhe main algorifhm for ASD of a random mafrix A for n = 20, where fhe number of bifs of 
precision used is sfandard Maflab single precision of 32 bifs, and require outpuf precision of 
5 = 10“^. We fhen run 50 • n ■ logu ~ 2500 copies of fhe fiber procedure Filter {A, m, 5). Then, 
we iferafe over a sequence of values for fhe infeger M, i.e. fhe number such fhaf such fhaf by 
Theorem 32 a uniformly random infeger m ^ U\M] separafes any eigenvalue of a perfurbed 
version of A w.h.p. 

One can see fhaf while fhe analyfic behavior requires M fo be af leasf here we see fhaf 
even M = already suffices for fhe algorifhm fo refum an approximafe eigenvector for each 
eigenvalue w.p. very close fo 1. We conjecfure, hence, fhaf implemenfing fhe algorifhm will 
resulf in far beffer pracfical run fime fhan fhe analytical bounds we provide here. 
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Figure 1: Algorithm Behavior 
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